
function MainStep1_ParallelCandi
%
% 并行计算LLO出发反向积分至LEO附近的近地点庞加莱图
%
% 输出：
% candi_LEO2LLO.mat
%
% 作者：张晨
% 邮箱：chenzhang@csu.ac.cn
% 单位：中国科学院空间应用工程与技术中心，空间探索室
% 时间：2022年02月24日
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear; close all; clear path;
% restoredefaultpath;
warning off;

% 添加模块
parentFolder = fileparts(pwd);
addpath(genpath([parentFolder , '\Basic']));
addpath(genpath([parentFolder , '\Bcr4bp']));
addpath(genpath([parentFolder , '\Cr3bp']));

% 载入参数
aux = Bcr4bp_Aux;
t_star = aux.dim.EMRot_t;

% --------------------------- 并行计算LLO - LEO近地点庞加莱图 ---------------------------
% 轨道数
calNum = 1e4;

% 结果预分配
resultMtx = cell(calNum , 1);

tic
WaitMessage = parfor_wait(calNum);
parfor iLoop = 1 : calNum
    % for iLoop = 1 : calNum

    WaitMessage.Send;

    % 随机生成LI（Lunar Injection）初值
    [tauf , xxf] = Bcr4bp_RandLI2D('emrot' , 1 , [0.65 , 0.65] , aux);

    % -------------------------------- 步骤1：近月点 -> WSB --------------------------------
    try
        options = odeset('RelTol' , aux.tol , 'AbsTol' , aux.tol , 'Events' , @Bcr4bp_EventEMRot2DEarthAp);
        [seg1_tt , seg1_xx , seg1_te , seg1_xe , seg1_ie] = ode113(@Bcr4bp_EqmEMRot2D , ...
            [tauf , tauf - 80 * 86400 / t_star] , xxf , options , aux);
    catch
        continue;
    end

    % 没有远地点，跳出
    if isempty(seg1_ie)
        continue;
    end

    % -------------------------------- 步骤2：WSB -> 近地点 --------------------------------
    try
        options = odeset('RelTol' , aux.tol , 'AbsTol' , aux.tol , 'Events' , @Bcr4bp_EventEMRot2DEarthRp);
        [seg2_tt , seg2_xx , seg2_te , seg2_xe , seg2_ie] = ode113(@Bcr4bp_EqmEMRot2D , ...
            [seg1_te , seg1_te - 80 * 86400 / t_star] , seg1_xe' , options , aux);
    catch
        continue;
    end

    % 没有近地点，跳出
    if isempty(seg2_ie)
        continue;
    end

    %     % 【测试，画EM轨道】
    %     Bcr4bp_ScnEMRot(seg1_tt(1) , aux);
    %     plotOrb2D(seg1_xx , 'b' , 2);
    %     plotOrb2D(seg2_xx , 'r' , 2);

    %     % 【画SB1轨道】
    %     [tt1_SB1Rot , xx1_SB1Rot] = Bcr4bp_EMRot2SB1RotMtx2D(seg1_tt , seg1_xx , aux);
    %     [tt2_SB1Rot , xx2_SB1Rot] = Bcr4bp_EMRot2SB1RotMtx2D(seg2_tt , seg2_xx , aux);
    %     Bcr4bp_ScnSB1Rot(tt1_SB1Rot(1) , aux);
    %     plotOrb2D(xx1_SB1Rot , 'b' , 2);
    %     plotOrb2D(xx2_SB1Rot , 'r' , 2);

    %     % 【测试2】
    %     options = odeset('RelTol' , intTol , 'AbsTol' , intTol);
    %     [seg_tt , seg_xx] = ode113(@Bcr4bp_EqmSB1Rot2D , [seg2_te , tau0] , seg2_xe' , options , aux);
    %     plotOrb2D(seg_xx , 'g' , 2);
    %     seg_xx(end , :) - x0

    % --------------------------- 保存转移轨道初值 ---------------------------
    tau0 = seg2_te;
    xx0 = seg2_xe;
    resultMtx{iLoop} = [tau0 , xx0 , tauf , xxf];

end
WaitMessage.Destroy;
toc

% ------------------------------ 保存结果 ---------------------------------
resultMtx(cellfun(@isempty , resultMtx)) = [];
candi_LEO2LLO = zeros(size(resultMtx , 1) , size(resultMtx{1} , 2));
for iLoop = 1 : size(resultMtx , 1)
    candi_LEO2LLO(iLoop , :) = resultMtx{iLoop};
end

save candi_LEO2LLO candi_LEO2LLO

fprintf('candi_LEO2LLO 计算完毕！ \n')
fprintf('轨道初值: %0.0f \n' , calNum)
fprintf('LLO反向积分至LEO轨道数: %0.0f \n' , size(candi_LEO2LLO , 1))

end
